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We develop a stochastic Gross-Pitaveskii theory suitable for the study of Bose-Einstein condensation in a 
rotating dilute Bose gas. The theory is used to model the dynamical and equilibrium properties of a rapidly ro- 
tating Bose gas quenched through the critical point for condensation, as in the experiment of Haljan et al. [Phys. 
Rev. Lett., 87, 21043 (2001)]. In contrast to stirring a vortex-free condensate, where topological constraints re- 
quire that vortices enter from the edge of the condensate, we find that phase defects in the initial non-condensed 
cloud are trapped en masse in the emerging condensate. Bose-stimulated condensate growth proceeds into a 
disordered vortex configuration. At sufficiently low temperature the vortices then order into a regular Abrikosov 
lattice in thermal equilibrium with the rotating cloud. We calculate the effect of thermal fluctuations on vortex 
ordering in the final gas at different temperatures, and find that the BEC transition is accompanied by lattice 
melting associated with diminishing long range correlations between vortices across the system. 

PACS numbers: 03.75.Kk,03.75.Lm,05.70.Fh,05.10.Gg, 



I. INTRODUCTION 

Rotating Bose-Einstein condensation was first realized by 
Haljan et al.at JILA |[lt], who grew large vortex lattices by 
evaporatively cooling a rotating thermal cloud. The evapo- 
ration is performed in a prolate trap, along the axis of rota- 
tion so that escaping atoms make large axial excursions but 
do not stray far from the rotation axis. Escaping atoms carry 
away significant energy but little angular momentum, increas- 
ing the angular momentum per trapped atom. Using this cool- 
ing mechanism, which requires typically ~ 50s of evaporative 
cooling and an extremely cylindrical trap, very large vortex 
lattices containing more than 300 vortices have been created 
i^]. This nucleation mechanism is distinct from more stud- 
ied vortex formation mechanisms H H Isllej ItHsI lol [Toj [TUl 
because, rather than manipulating a condensate to generate 
vortices, vortices can be formed at the onset of condensa- 
tion as the rotating gas reaches the critical temperature (Tc) 
for Bose-Einstein condensation. This mechanism thus has 
connections with the Kibble-Zurek lfl2l [isll and Berezinskii- 
Kosterlitz-Thouless 
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mechanisms, as it is funda- 
mentally caused by spontaneous symmetry breaking and ther- 
mal fluctuations near a phase transition. 

In this work we study the dynamics of vortex lattice for- 
mation during rotating evaporative cool ing using a stochastic 
Gross-Pitaevskii equation (SGPE) HI [11, [H]. Our central 
aim is to include rotation, thermal fluctuations, and two-body 
collisions, by separating the system into a condensate band 
coupled to a (possibly rotating) thermal reservoir. The con- 
densate band is described by a generalized mean field theory 
known as the truncated Wigner method ll20l l2ll I22I1 that in- 
cludes the projected classical field method ll23l |24 I25I1 as 
a special case in the regime of high temperatures, while the 
non-condensate band formalism is based on quantum kinetic 
theory ll26ll . 

We also provide expressions for the dissipation rates of the 



theory under the condition that the high-energy component 
is approximately Bose-Einstein distributed, and details of our 
numerical algorithm that is central to the utility of the SGPE 
approach. These are the first SGPE simulations of vortex for- 
mation in a trapped Bose gas. 

II. SYSTEM AND EQUATIONS OF MOTION 

The trapped system is divided into a condensate band of 
states with energy beneath the cut-off and the remain- 
ing non-condensate band of high-energy states. The formal 
separation of the Hamiltonian into bands according to energy 
allows different methods to applied in the two bands. The 
condensate band is treated using the truncated Wigner phase 
space method ll20ll2l[ | which is valid when the modes are sig- 
nificantly occupied ll27ll . Within this approach the conden- 
sate band includes the condensate and the most highly occu- 
pied excitations. The non-condensate band is comprised of 
many weakly occupied modes and is treated as approximately 
thermalised; the effect of this band is to generate dissipative 
terms in the equation of motion for the condensate band. The 
high-energy states play the role of a grand canonical thermal 
reservoir for the low energy condensate band. In this section 
we provide a minimal overview of the formalism and simply 
state the final equations of motion to be used in this work. 
Refs. ll^lM 

contain the essential formalism, and Ref. Il28ll 
gives necessary modifications for the rotation system, along 
with a number of corrections to the derivation of the equations 
of motion. 



A. Hamiltonian and separation of the rotating system 

The cold-collision regime for bosonic atoms of mass m 
is characterized by the interaction parameter u = Antra I m. 
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FIG. 2: Interactions between the condensate and non-condensate 
bands. In (a) two non-condensate band atoms collide. The collision 
energy is transferred to one of the atoms with the other passing into 
the condensate band. In (b) a condensate band atom collides with a 
non-condensate band atom with no change in condensate band pop- 
ulation. The time reversed processes also occur. 



FIG. 1 : Separation of the system into condensate and non-condensate 
bands. States beneath the energy E^, defined in the rotating frame, 
form the condensate band, while states above Er form the non- 
condensate band consisting of high-energy, approximately ther- 
malised atoms. In the s-wave scattering regime high-energy states 



{Es 



) are also eliminated from the theory IT^ . 



where a is the s-wave scattering length II29I1 . The second- 
quantized many-body Hamiltonian for a dilute Bose gas con- 
fined by a trapping potential V{x, t) in the rotating frame de- 
fined by angular frequency £2 = Qz is i/ = i/sp + Hj, where 



d\ iA"'"(x) -— - ni, + v(x, f)U(x), 



and 



u C 
"'-2 J 



(1) 



(2) 



The orthogonal projectors that separate the system are defined 
with respect to the single particle basis of the trap. This pro- 
vides a convenient basis because (a) The many body spectrum 
becomes well approximated at high-energy by a single par- 
ticle spectrum, allowing the cutoff to be imposed in this ba- 
sis 1I24I1 . and (b) the numerical method used to solve the final 
equations of motion uses the single particle basis Isoll . 
The condensate band projector takes the form 



^ = |n><n|, 



(3) 



where the bar denotes the high-energy cut-off and n represents 
all quantum numbers required to index the eigenstates of the 
single particle Hamiltonian. To minimize notation we will use 
!P both for the operator and for its action on wave functions; 
the meaning will always be clear from the context. The field 
operator is decomposed as 



iA(x) = 0(x) + -A^c(x) = H{x) + Q4>(x), 



(4) 



where the non-condensate field {ff^cOO = <3iA(x) describes the 
high-energy thermal modes, and the condensate band is de- 
scribed by the field ^(x) = !P0'(x). In the spatial representation 
the eigensolutions of H^pYa = en^^n, denoted Ya{x) = (x|n). 



generate the action of the projector on an arbitrary wave func- 
tion *P(x) 

P^'ix) = Yj ^"(x) r Y»nz)- (5) 

n ^ 

The condensate band field theory generated by such a de- 
composition has a nonlocal commutator 



[0(x),0"'"(y)] =5c(x,y) 



(6) 



where 6c(x,y) - {x|!P|y). In the harmonic oscillator repre- 
sentation of a one dimensional system, the approximate delta 
function has a position dependent width, being narrowest in 
the center of the trap. This reflects the fact that trap wave- 
functions become smoother near their semi-classical turning 
points. Note that since !PP = P, !P0(x) - 0(x); put another 
way, this means that (5c(x, y) is a true Dirac-delta function for 
any condensate band field. 

There are two physical processes arising from the inter- 
actions between condensate and non-condensate bands, re- 
ferred to as growth and scattering UtI IT9I l28ll . and shown 
schematically in Figure [2] The growth terms correspond to 
two-body collisions between condensate and non-condensate 
atoms whereby particle transfer can take place. The scattering 
describes inter-band collisions that conserve the populations 
in each band but involve a transfer of energy and momentum 
between bands. The non-condensate band is described semi- 
classically in terms of the Wigner function 

F(x, K) ^ J d\ <(A]vc (x + I) he (x - e'^ \ (7) 

In this work the trapping potential assumes the cylindrical 
form 



V{x,t)=-{ajy+ulz\ 



(8) 



where and lo- are the radial and axial trapping frequencies 
respectively. The non-condensate band Wigner function for a 
Bose gas confined in a cylindrically symmetric trap in equi- 
librium in the rotating frame at frequency Q, with chemical 
potential /i, is 



F{x, K) : 



1 



exp [(^w(x, K) - ^i)|kBT] - 1 ' 



(9) 
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where 



nio(x, K) =— ■ (x X K) + V(x), 

2m 

= ^(K-maxx/^)2 + yeff(x), 
Zm 



(10) 



defines the effective potential as VeffCx) - m{ojf. - Q^)^r^/2 + 
m(x)?z^/2, and r - (x^ + y^)'^^ is the distance from the sym- 
metry axis of the trap. Since = K - mil xx/his the rotat- 
ing frame momentum, we can work in terms of this variable, 
defining 

hcaix,K) = ;!w,(x,K,.) = — — ^ + yeff(x). (11) 
2m 

Hereafter we drop the subscript for brevity since we al- 
ways describe the system in this frame. In the semiclas- 
sical approximation f(u, K) is only significant in the non- 
condensate region Rnc^ which is the region of phase space 
where ^(jj(u, K) > £«. 



B. Stochastic Gross-Pitaevskii equation 

Following the theory in Refs. lfT7i[T9ll28ll . the equation of 
motion for the density operator describing the total system 



(12) 



is reduced to a master equation for the condensate band den- 
sity operator, which is the trace of p over the non-condensate 
band degrees of freedom pc = ir^c ip)- The master equation 
is then mapped to a Fokker-Planck equation for the Wigner 
distribution describing the condensate band. In the truncated 
Wigner approximation, an equivalent stochastic equation of 
motion can be obtained that allows the dynamics to be effi- 
ciently simulated numerically jsill . Formally, there is a cor- 
respondence between symmetrically ordered moments of the 
field operator 0(x, t) and ensemble averages over many tra- 
jectories of the projected field a{x, t) = y,„an{t)Y„{x) that 
evolves according to the SGPE to be defined below ll28ll . For 
example, for the density 

{a\x)a{x))w = i<0'"(x)0(x) + 0(x)0'"(x)), (13) 

where ( )w denotes stochastic averaging over trajectories of 
the SGPE. 

To state the SGPE we require some expressions from mean 
field theory, in particular the GPE obeyed by a{x) in the ab- 
sence of any thermal atoms in the non-condensate band. The 
mean field Hamiltonian for a{x) corresponding to H is 



2m 



-D.L- + V{x)\a{x) 



+ \^d\ \a{x)\\ 
with atom number given by A^gp — J ^''x |a(x)|^. 



(14) 



The growth terms involve the Gross-Pitaevskii equation 
found in the usual way as 



^^(x) 6Hcp , , 

in — - — = PLcpaix) 



dt fo*(x) 



2y72 



2m 



- Q.L, + V{x) + u\a{xT a{x) 



(15) 



that defines the evolution operator Lgp- To obtain this equa- 
tion of motion we have used the projected functional calcu- 
lus ifigjl that plays the same role here as ordinary functional 
calculus in non-projected field theory. 

The scattering process couples to the divergence of the con- 
densate band current 



jGp(x) = A([Va*(x)]a(x) - a*(x)Va(x)) 



2m ^ 

-{nxx)\a(x)\\ 



(16) 



where the second line is a rigid body rotation term arising 
from the transformation to the rotating frame. The limit 
jGp(x) - gives the laboratory frame velocity field v = Q x x, 
corresponding to the irrotational system mimicking rigid body 
rotation. Persistent currents around the vortex cores are re- 
sponsible for maintaining irrotationality. 



1. Full non-local SGPE 

We finally arrive at the Stratonovich stochastic differential 
equation (SDE) IH 



{S)da(x, f) = - -PLGpa(x)dt 
n 

G(x), 



+ r d'x' M(x - x')V ■ jGp(x')flff 



(17) 



+ ia(x)dWM{x, t) 



where the growth noise is complex, and the scattering noise is 
real. The noises are Gaussian, independent of each other, and 
have the non-zero correlations 

{dWl{x,t)dWc{x',t)) = 2G{x)6c{x,x')dt, (18) 
{dWM{x,t)dWM{x' J)) = 2M(x-x')^ff. (19) 

The growth and scattering amplitudes, G(x) and M(x), are cal- 
culated in AppendixlAlfor the case where the non-condensate 
band is described by a thermal Bose-Einstein distribution. 

The two validity conditions for this equation are (i) high 
temperature: ho) <K kgT where a> is the characteristic system 
frequency, and (ii) significant occupation of modes in the con- 
densate band, necessitated by the truncated Wigner method. 
The first fine of Eq. [17] is the PGPE developed by Davis et 
al ll24ll that, in essence, is the GPE formally confined to a 
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low energy subspace of the condensate band. The conden- 
sate growth terms of the second Hne are closely connected to 
phenomenological dissipative theories of BECs based on the 
GPE ll32l [33I I34I1 but here have a specific form given in Ap- 
pendix|A] The remaining terms describe scattering processes 
that transfer energy and momentum between the two bands; 
these terms where first derived in Ref. |fl9ll. and have not ap- 
peared in previous SGPE theories. 

Irrespective of the form of G(x) and M(x), the equilibrium 
state for the Wigner distribution of the system has the grand 
canonical form Iil9il 
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Ws exp 



^A^GP ~ ^GP 



(20) 



that is generated by the stochastic evolution of Eq. (fTTT i. So- 
lutions of the SGPE also known to satisfy a set of general- 
ized Ehrenfest relations which in the steady state reduce to a 
statement of the fluctuation-dissipation theorem for the grand 
canonical Bose gas llssll . 



FIG. 3; (a) Variation of ideal gas BEC transition temperature with an- 
gular frequency of rotation fi, for N = 1.3 x 10^ atoms, (b) Angular 
momentum per particle (solid line), for the same number of atoms at 
temperature T = Tc, with {L-)l{TiN) ± (r(L-)/1i (dashed lines) repre- 
senting the size of fluctuations. Both the angular momentum per par- 
ticle and its fluctuations diverge at at the critical rotation frequency. 



2. Simple growth SGPE 



ni. SIMULATIONS OF ROTATING BOSE-EINSTEIN 
CONDENSATION 



A simpler equation that reaches the same equilibrium state 
can be obtained by noting that the scattering terms in Eq. ( fTTl i 
do not directly change the condensate band population. Such 
terms have been included in a kinetic theory of Bose-Einstein 
condensation ll36ll . where they were found to play a role in the 
early stages of condensate growth by causing atoms in many 
modestly occupied states to merge more rapidly into a single 
highly occupied state which then dominates the growth pro- 
cess. Within the SGPE theory, the scattering terms generate 
an effective potential from density fluctuations. This can be 
seen by noting that upon neglecting the growth and scattering 
terms, that are relatively small compared to the mean field 
evolution, the continuity equation for the condensate band 
gives V ■ jGp(x) = -c?«Gp(x)/(9f. Consequently, since M(x) is 
purely real, the deterministic part of the scattering terms gen- 
erate an effective potential that will tend to smooth out density 
fluctuations. 

Upon neglecting these terms the resulting stochastic equa- 
tion is much easier to integrate numerically as it does not re- 
quire the non-local noise correlations of the full equations. 
The simulations in this work are performed by numerically 
integrating the simple growth SDE 



da(x) 'PLcpaix)dt 

ft 



^{-^{fi- Lcp)aix)dt + dWyix, t) 



(21) 



using the mixed quadrature pseudo-spectral method detailed 
in Appendix IB] The rate y is given by Eq. ( lAl lb . and the 
noise is complex Gaussian with non-vanishing correlator 

{dW;ix, t)dWyix', 0) = 2y6cix, x)dt. (22) 



The central question of interest for this work regards the 
formation of vortices during condensation. In particular, does 
a vortex-free condensate form and then vortices enter it from 
the edge, or does condensation proceed into a state with vor- 
tices already present? That this is non-trivial to answer is eas- 
ily seen from the single particle energy spectrum of the cylin- 
drically symmetric harmonic trap in three dimensions 



\) - nQ.1 + Tioj^im + 111), 



(23) 



where n, /, m are the radial, angular and axial quantum num- 
bers. In the absence of interactions BEC will form in the 
ground state mode with energy eooo - fi<^r + ficozl^- Vortices 
arise from occupation of states with nonzero angular momen- 
tum and the positive angular momentum part of the spectrum 
behaves as h{LOr — Q)/ leading to near-degeneracy of positive 
angular momentum modes as Q — > w, . States with angular 
momentum {L,} - M > sue increasingly easy to excite as 
the rate of rotation increases and in the interacting gas there is 
a trade-off between the interaction and rotational contributions 
to the energy due to the presence of a vortex that is known to 
favor the creation of vortices in the condensate once the sys- 
tem is rotating faster than the critical angular frequency, Q.c, 
for energetic stability of a vortex iIstIi . For a rapidly rotating 
BEC transition (O » Qc), we can expect that vortices will 
play a dominant role at all stages in the BEC growth process. 
It is then possible that atoms condense into a vortex-filled 
state. In general the vortex configuration can be disordered. 
It is not required, for example, that the vortices minimize the 
system energy by forming a regular Abrikosov lattice at all 
stages of condensate growth. 
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FIG. 4: (color online) The Wigner function for the rotating ideal Bose gas used as the initial state for our simulations of rotating Bose-Einstein 
condensation. Parameters are Tq = 12nK, SIq = 0.979a;,-, = O.Seooo- (a) Mean density {(^^(x)ij(x)> computed by averaging over 5000 samples 
of the Wigner function. A single sample is shown in (b) with vortices (+ have positive angular momentum) shown for r,, < 3ro. (c) Condensate 
wave-function according to the Penrose-Onsager criterion. The one body density matrix used here is constructed from 500 samples. Due to 
the near-degeneracy of angular momentum modes averaging must be taken over many samples to completely remove all vortices from the 
proto-condensate. In this example the proto-condensate contains ~ 50 atoms and many vortices are visible. In all plots the colormap is scaled 
to represent the entire range of density in the data. 
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FIG. 5: Condensate number A'o during a quench from [}!(,, Tg] = 
[0.5eooo, 12nK] to [p, T] = [3.5eooo, InK] at fixed rotation frequency 
f2o = 0.979a;, . The condensate number is found by diagonalising the 
one-body density matrix found from short time averaging single runs 
of the SGPE (see text). The final condensate number is suppressed 
at higher temperature. 



A. Properties of the rotating ideal Bose gas and the 
non-condensate band 

An important effect of rotation in Bose-Einstein conden- 
sation is the reduction of Tc llssi [39tl . that, along with the 
slow one-dimensional evaporative cooling, accounts for the 
long evaporation time in the JILA experiments. The transition 
temperature for the ideal Bose gas is modified by the effective 
radial frequency in the rotating frame, cl>j_ = ^co^ - DP-, so 
that 



2\l/3 



(24) 



where r° = 0.94-h(Liyj.Lij,y^^N^^^ /kg is the nonrotating transi- 
tion temperature for trapped atoms. 

In this work we model the non-condensate region assuming 
the high-energy atoms are quickly thermalised into a Bose- 
Einstein distribution, that we approximate as non-interacting 
since the density is usually small at high energies. We can 
find analytical expressions for the properties of the noninter- 
acting gas, including the effect of the cutoff, by defining an 
incomplete Bose-Einstein function 



gyiz,y)=— dxx'-'Y^ze-'')', 

r(v) Jy 



y z' r(v,yD 

r(v) ' 



(25) 



where r(v, x) = dy y^'^ e^-' is the incomplete Gamma func- 
tion. In analogy with the reduction to an ordinary Gamma 
function r(v,0) = r(v), we have gy(z,0) = gy(z) = S^i z'/T, 
reducing to the ordinary Bose-Einstein function. We then find 
for the non-condensate band density, for example 



= ^ Ve^^-^-W] dkk^e-"' 



(26) 



where fi^Kii(x)^/2m = max{ER - V'efF(x), 0). Changing inte- 
gration variable toy = n/3trk^/2m gives 



riNC 



(X) = gy2 (/^-''-«] ,;8/z'/:«(x)V2m]) /Al^, (27) 



where A^g - -sjTjr^Jmk^ is the thermal de Broglie wave 
length. Setting Er = 0, we recover the standard form for the 
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semi-classical particle density of the ideal gas 112911 . The total 
number in the non-condensate band is given by 



(28) 



where cj - (w-w^)'^^ is the geometric mean frequency in the 
rotating fame. In this way the usual semi-classical expressions 
can be generalized to include a cutoff in terms of the incom- 
plete Bose-Einstein function. 

The angular momentum and its fluctuations 
can also be calculated. We now write (A) = 
{2nY^ ^ (fix ^ (P^ A(x, k)F(x, k), and the variance per 

particle as (T(Af = {A^)INnc - {{A)INNcf. In ro- 
tating frame coordinates the angular momentum is 
L~ = fi(xky - yk^) + mQ.{x^ + y^), which in equilibrium 
gives the rigid body relation between angular momentum and 
rotational inertia (L-) = mQ.{x- + y^). From this expression 
we find the angular momentum is 



(L.) = in 



(29) 



where w' = (w^w^)'^''. 
momentum is 



Similarly, the mean squared angular 



2^^g5(g^^yg£«) 



1 + 



4Q2 



(30) 



where a) = {a)^a>\)^^^. In Figure[3]we give an example of the 
the variation of Tc, and the angular momentum per particle 
together with its fluctuations as a function of the angular fre- 
quency of rotation. The suppression of Tc by the rotation is 
seen to be associated with diverging angular momentum and 
angular momentum fluctuations near the critical frequency for 
stable rotation in the harmonic trap. 



B. Numerical procedure 

The long cooling time and the large changes in length scales 
in the axial and radial dimensions during cooling makes mod- 
eling the full condensation dynamics challenging. However, 
as the gas cools and spins up, the system becomes progres- 
sively more two-dimensional. We model the system by con- 
sidering the two dimensional projected subspace of the quasi- 
two dimensional system. We thus use the SGPE of Eq. |2T| 
with rescaled interaction parameter U2d = AntP'a/mL-, where 
in this work we take the thickness through the oblate system 
as L- mil)-, obtained by projecting onto the ground 

state of the z-dimension. Note that a consequence of the ap- 
proximation G(x) ^ y is that the growth rate is unchanged by 
reduction to the two dimensional effective equation. 

We take as our starting point the gas after it has been evap- 
oratively cooled to a rapidly rotating state above the transition 
temperature. We then simulate the dynamics of a rapid quench 
below threshold. The procedure is summarized as follows: 

(i) Initial states are generated by sampling the grand 
canonical ensemble for the ideal Bose gas for initial pa- 
rameters Tq, D.Q and yUo- 



(ii) To approximate a rapid quench, the temperature and 
chemical potential of the non-condensate band are al- 
tered to the new values T < To and yU > /io, while the 
rotation of the cloud Q is held fixed at Qq, consistent 
with axial evaporation. The condensate band field is 
then evolved according to the SGPE. 

Specifically, we aim to model a possible run of the JILA 
experiment consisting of a system of '^^Rb atoms held in a 
magnetic trap with frequencies (w.,(y,.) = 27r(5.3, 8.3)Hz, 
initially at temperature To = 12 nK with chemical poten- 
tial //o = O.Seooo, and rotating with angular frequency = 
0.979wr. We examine a range of final quench temperatures 
T = [1, 2, 3,4, 5,7, 9, 11] nK with constant chemical potential 

- 3.5eooo- We note the ratio Q/jj = 0.28 <*: 1, where a 
value approaching unity indicates the onset of the mean field 
quantum Hall regime. The ratio to- / ^jaP^ - D?- ~ 3.13 » 1, 
indicates the system is oblate but not highly two-dimensional. 
We have found that simulations for parameters that are highly 
oblate are very challenging numerically, so we consider as 
our first application of the method these more modest param- 
eters. For the results presented here we use a cut-off value 
of En = SficiJr, giving a condensate band containing 291 sin- 
gle particle states. To check the cut-off independence of our 
results we have also carried out simulations for Er = 6^Wr 
(not shown here). We find qualitative agreement for single 
trajectory dynamics and quantitative agreement (at the level 
of a few percent) for the ensemble averaged condensate frac- 
tion shown in Figure[TO] We are thus confident that our results 
are cut-off independent. For computational reasons we choose 
constant dimensionless damping hy/keT - 0.01 for all simu- 
lations. Although we have derived a form for the damping rate 
-y given in AppendixlAl the values of y arising from Eq. (lAl II) 
for our parameters are orders of magnitude smaller then our 
chosen value, so we have chosen a fixed damping such that 
the growth rate is comparable to the JILA experimental time 
of 50s. This affects the specific time-scale for growth, but 
not the qualitative features of the simulations which merely 
require that the damping is much slower than any other time- 
scale of the system. We emphasize that the equilibrium prop- 
erties are independent of y. 

Since our numerical simulation method employs a pro- 
jected eigenbasis of the associated linear Schrodinger prob- 
lem, the initial Bose-Einstein distribution is easily sampled. 
The grand canonical Wigner distribution for a non-interacting 
trapped Bose gas distributed over a set of single particle states 
with frequencies w,,/ is 



\a„iy 



N„i + 1/2, 



(31) 



where a{x) - 2„;a„/F„;(x), the modes are defined in Eq. 
and 



1 



,{tioj„,-no)lkBTa _ \ ' 



(32) 



The energy cut-off is denoted by the upper limit of the product 
notation. The single particle spectrum of the condensate band 
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FIG. 6: (color online) Condensate band density for a single trajectory of the SGPE. (a) Initial state for 1.3 x 10' ^'^Rb atoms at To = 12nK, 
with ^0 = 0.5eooo> and ^^o = 0.979w,-. At time t = the non-condensate band is quenched to T = InK, fi = 3.5eooo> and fl = Slg. (b)-(d) The 
condensate band undergoes rapid growth, (e)-(f) At this low temperature the vortices assemble into a regular Abrikosov lattice. 



is restricted by the cutoff to 

nuj„i = hur{2n + \l\ - QZ + 1) < Er. (33) 

The distribution is sampled as a„i = '\/N„i + 1/2 ri„i, where 
T]„i are complex Gaussian random variates with moments 
{rinirjmq) = W„,J]mq) = 0, and {T]l,J],„q) = 5„,„5/,. We simulate 
the time dynamics of the quenched system using the numeri- 
cal algorithm described in Appendix 151 

We note here that application of the Penrose-Onsager cri- 
terion for obtaining the BEC ll30l l40ll is hindered in this sys- 
tem because each run of the SGPE theory leads to vortices 
in different locations, and ensemble averaging the trajectories 
generates spurious results for {i^' (x)0'(x')). Specifically, the 
condensate number is significantly reduced below the conden- 
sate number that we find using a different technique of short- 
time averaging single trajectories. Generalizing the Penrose- 
Onsager criterion for BEC to multi-mode mixed states is not 
pursued in this work. Here we have taken the view that a given 
trajectory samples a sub-ensemble of the full range of quan- 
tum evolution, consistent with spontaneously broken symme- 
try in the initial conditions. We present single trajectory re- 
sults of our simulations of the SGPE that give a qualitative 
indication of what could be seen in individual experimental 
measurements of the atomic density time series. 



C. Initial states and dynamics 

In Figure|4]we show the density of the condensate band for 
the thermal gas above the BEC transition. Individual samples 
of the Wigner function are seen to contain phase singularities. 
We also construct and diagonalize the one body density ma- 
trix, according the the Penrose-Onsager criterion for BEC. By 
averaging 500 samples from the Wigner function we find that 
the largest eigenvalue is A^o ~ 50 and the associated eigen- 
vector contains several phase singularities. Using 5000 sam- 
ples washes out the singularities and gives the noninteracting 
ground state; the point here is that the high degree of degen- 
eracy means that we must average over very many samples 
to remove the singularities. Put another way, this means that 
individual trajectories always contain many vortices, and the 
motion of these vortices destroys long range coherence asso- 
ciated with BEC. 

In Figure |5] we plot the condensate number as a function 
of time after the quench for a range of different quench tem- 
peratures. The condensate is calculated by constructing the 
one body density matrix for the sub-ensemble via short time 
averaging individual trajectories ifsoll . We have found that by 
averaging over 50 samples taken uniformly over an interval of 
2.5 radial trap periods, we obtain sufficiently good statistics 
to determine the condensate number A^o, given by the largest 



8 




FIG. 7: (color online) Condensate mode calculated from the Penrose-Onsager criterion via short time averaging, for the simulation data of 
Fig|6] (a) Initial state, (b)-(d) After the quench the condensate undergoes a rapid growth into a highly disordered phase, (e)-(f) The system 
undergoes a prolonged stage of ordering into a periodic vortex lattice. At this temperature the condensate and condensate band are almost 
identical in the final equilibrium state. 



eigenvalue of (i^' (x)i/'(x')). The main point to note from this 
result is that at higher temperatures the absolute condensate 
number is slightly suppressed. However, as we will see in the 
next section, the condensate /racf/on is strongly suppressed, 
due to the increasing atom number in the non-condensate band 
for higher temperatures. 

A typical time series is shown in Fig.|6] for the lowest tem- 
perature considered in Fig. |5] T = 1 nK. The initial sample 
contains many phase singularities, including a few that are 
circulating counter to the cloud rotation. Counter-circulating 
vortices only tend to arise in the initial state, are thermody- 
namically very energetic, and are quickly annihilated by pos- 
itive vortices during post-quench dynamics. The condensate, 
shown in Fig. |7l forms with many vortices already present, 
but in a highly disordered configuration. The vortices then 
assemble into a regular triangular lattice. 

At low temperatures the equilibrium state of the system 
consists of a an Abrikosov vortex lattice with small ther- 
mal fluctuations of the vortices about their average positions. 
The time dynamics for the highest temperature of Fig. |5] 
r = 11 nK, are shown in Fig. [8] We see that growth pro- 
ceeds into a disordered state for which periodic order remains 
frustrated by significant thermal fluctuations. 



D. Equilibrium states 

We now consider equilibrium properties of the rotating sys- 
tem, in the grand canonical ensemble generated by evolu- 
tion according to the SGPE. In Fig. [10] we plot the conden- 
sate fraction of the final state for the temperatures used in 
Fig. El The total atom number in the system Nt is calculated 
as Nt = Nc + Nnc, where N^c is given by Eq. ( |28] ) evaluated 
for our chosen energy cutoff, and Nc is the mean total number 
of atoms in the condensate band. 

We note that the ideal gas relationship between condensate 
fraction and temperature, No/N = 1 - (T/Tc)^, is unaltered 
by rotation, provided the rotating frame transition tempera- 
ture, given by Eq. [24l is used in place of the non-rotating 
transition temperature. The rotating ideal gas curve is plot- 
ted for comparison, and we see that the condensate fraction 
predicted by the SGPE simulations is suppressed relative to 
the ideal gas result. Although we are unable to simulate the 
nonrotating system with our two-dimensional algorithm, since 
the ideal gas behavior is universal with respect to rotation we 
also compare with the experimental data for the non-rotating 
Bose gas ll4lll . Our simulations show a condensate fraction 
significantly suppressed relative to the experimental data for 
the non-rotating system. 
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FIG. 8: (color online) Condensate band density for a quenched rotating Bose gas. (a) Initial state for 1.3 x 10 Rb atoms at To = 12nK, with 
fio = 0.56000, and CI = 0.979a;,-. At time / = the non-condensate band is quenched lo T = llnK, ft = 3.5eooo> and CI = Clg. (b)-(d) The 
condensate band undergoes rapid growth, (e)-(f) After growth, vortex order is suppressed by thermal fluctuations, frustrating the formation of 
a periodic vortex lattice seen at the much lower temperature of Fig. [6] 



This difference presumably comes about from two effects. 
The most obvious is the loss of long range order caused 
by thermal motion of vortices that are absent from the non- 
rotating system. A less obvious and possibly less significant 
effect is that the reduction of condensate number by interac- 
tions is more pronounced for a rotating gas, due to the exclu- 
sion of atoms from vortex cores. Our highest temperature sim- 
ulation, at r = 11 nK, is right on the edge of Bose-Einstein 
condensation, and this is reflected in the lack of long range 
order seen in Fig. [8] 

The loss of long range order can be seen more clearly by 
looking at the relative positions of vortices in the condensate 
wave-function at different temperatures 14211 . In Figure [TTIwe 
plot the condensate wave function, along with the histogram 
of the separation of each pair of vortices in the system. We 
only consider vortices with radius r,, < lOro. For the lowest 
temperature (T - I nK) the histogram is indistinguishable 
from the ground state histogram (not shown). The presence 
of many peaks for large separations reflects the long range 
spatial order in the system. As temperature increases, order is 
reduced starting at large vortex separation. At the highest tem- 
perature (r = 1 1 nK) we see that only two histogram peaks 
clearly persist, corresponding to nearest and next-to nearest 
neighbor order being preserved. 



IV. CONCLUSIONS 

We have developed and applied the SGPE theory to the ro- 
tating Bose gas. The formalism presented here can be equally 
applied to the non-rotating case and in that limit extends the 
prior SGPE theory |[l3,[il,|2l by providing expressions for 
the dissipation rates. We have developed an efficient algorithm 
for numerically integrating the SGPE in a rotating reference 
frame and applied it to the case of quenching a rotating Bose 
gas below threshold, modeling the rapid-quench limit of the 
experiment at JILA [Ql] . 

Our simulations indicate that many vortices co-rotating 
with the thermal cloud are trapped in the emerging conden- 
sate, reflecting the large angular momentum of the system. In 
our calculations these primordial vortices are nucleated from 
vortices in the individual samples of the initial Wigner dis- 
tribution and from fluctuations arising from the condensate 
growth terms in the SGPE. The process can be regarded as 
stimulated vortex production, in contrast with the non-rotating 
case where vortices can appear spontaneously but the total 
vorticity in the BEC must be zero on the average ll43ll . 

For low temperatures the vortices order into a regular 
Abrikosov lattice as the system approaches thermal equilib- 
rium with the quenched thermal cloud. At sufficiently high 
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FIG. 9: (color online) Condensate mode calculated from the Penrose-Onsager criterion via short time averaging, for the simulation data of 
FiglD (a) Initial state, (b)-(d) After the quench the condensate growth into a highly disordered vortex configuration, (e)-(f) Suppression vortex 
order is reflected in the disordered condensate wave-function. At this higher temperature the condensate and condensate band are distinct, 
reflecting the disorder and the increased thermal fraction in the condensate band. 



temperatures lattice crystallization is frustrated by thermal 
fluctuations of the positions of vortex cores. This disorder 
is reflected in a loss of long range correlations in the relative 
positions of vortices that are most evident as the system ap- 
proaches Tc- 

We note that the concept of rotation has not been in- 
vestigated in the context of the Kibble-Zurek mechanism, 
which explains vortex nucleation during the phase transition 
in terms of spontaneous symmetry breaking and critical slow- 
ing down 14411 . Future quantitative studies of dynamics at the 
transition will allow for a deeper understanding of the conden- 
sation process, in particular the role of spontaneous symmetry 
breaking as a vortex nucleation mechanism in rotating Bose- 
Einstein condensation. 
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APPENDIX A: CALCULATION OF DISSIPATION RATES 

Reservoir interaction rates closely connected to those re- 
quired here have been previously computed within the context 
of the quantum kinetic theory (QKT) of BEC ll26ll . We now 
derive the form of the growth and scattering rates including 
the effects of the cutoff, the external trapping potential and the 
rotation of the thermal cloud. 



1. Growth 

a. General form 
From Eq. (84) of Ref. lfl9ll . the growth amplitude is 

G(x) = ^ £/^vG^"'(x,v,0) = Gi(x)H-G2(x) (Al) 

where 

^i(^^^7t4f r r r d'^yd'Y^^d^Y., 

{2nyn- J J Jr^^ (A2) 
xF(x,Ki)F(x,K2)Ai23(0,0), 
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FIG. 10: Equilibrium condensate fraction versus temperature. The 
boxes show the result of SGPE simulations for the rotating Bose gas 
at constant rotation frequency Q. = 0.979aJr. The solid line is the 
ideal gas expression, which is invariant with respect to Q. For com- 
parison the open circles are the exp erimental data for the non-rotating 
system reported in Reference 14111 . 



X F{x, Ki)F(x, K2)F(x, K3)Ai23(0, 0), 



(A3) 



and Ai23(k, e) = 6(Ki + K2 - K3 - k)5(wi + cd2 - (^3 - e/fi) 
conserves energy and momentum. Note that effect of working 
in the rotating frame is only to introduce the effective potential 
term -mO^r^ jl into the cut-off and Wigner functions in G(x), 
which is otherwise unaltered. We will take the reservoir to be 
in equilibrium with Bose-Einstein distribution at chemical po- 
tential // and temperature T . Integration over the momentum 
<5-function reduces Gi(x) to 



Gi(x) 



Sottt^ 



Vdr(x))(P+9) 



p.q=l 

K2dK2 e-0'^r+«^iWr/2m 

j dz6{z-mV,six)/h^KiK2), 



/-*CO 



K,dK 



(A4) 



where Eq. ( fTTT ) gives the expression for the momentum cutoff 
Kr(x) = V2m[£fl- Veff(x)]/S2 (A5) 

Changing to variables s = p/SfP'K^/lm, t = q/Sti^K^ /2m, and 
taking account of the delta-function, we obtain 



Gi(x) 



p,q=i 



pq 



X 



I ds e ' I dt e ', 



(A6) 



where s„„„(x) = p/^[Er - Veff(x)], and f„„„(x, i) = 
mml^q/3[ER - Ves(x)],(J3Veftix)fpq/4s^. We now show that 
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FIG. 1 1 : (color online) Condensate density and vortex position his- 
tograms. Top: T = \ nK, middle: T = 5 nK, bottom: T = 11 nK. 
The highly ordered Abrikosov lattice has its order degraded by ther- 
mal fluctuations for increasing temperature. 



the rate is composed of two regions: an inner region where 
the rate is spatially invariant, and an outer limit where there is 
a weak spatial dependence. 



b. Inner region 

From the from of the lower limits we find that wherever 
yeff(x) < 2Er/3, t„,i„{x,s) = i„„„(x) = qP[ER - Vsd^)]. 
Proceeding similarly for G2(x) we find the same condition to 
hold, so that in the inner region the growth rates are spatially 
invariant, taking the final form 



G'l" 
Gif 



4m 
Am 



(A7) 



(fl/tBr)2[ln(l -e^*''-^"')] , 
X 2 e'-/*<^-2£«) (^^[(J^(^'-'^'<\ l,r+ l]f , (A8) 



where the Lerch transcendent is defined as <b[z, s, a] = 
YjT=(i^^ l^'^ + ky . These terms are superficially similar to the 
condensate growth rates previously obtained using the quan- 
tum kinetic theory ll45ll . although the form of the cutoff de- 
pendence is different in the SGPE theory. Note that this form 
is independent of any collisional contribution to the effective 
potential. 
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c. Outer region 

In the region where 2Er/3 < Vssix) < Er the rates depend 
on position, albeit quite weakly. At the edge of the condensate 
band where VeffCx) = E^, we find 



p,q=\ 



yfpq 



-/TiOSVm) (A9) 



„«/J-£«)(p+9)+;3(/J-2£fi)r 



XK,{J3^l{p + r){q + r)\ 



(MO) 



where is the modified Bessel function of the second 

kind. 

In Fig.[T2]we compute the variation of Gi(x) with radius 
for a spherical parabolic potential. For illustration we have 
also neglected the collisional contribution to the effective po- 
tential. To show the maximum variation of the total rate we 
also plot (inset) the ratio Rc = (Gf" + G°"')/(Gf + G^') ver- 
sus nlksT . For all data we use Er = 3// which is typically 
a good approximation for the energy at which the spectrum 
of the many body system becomes single particle-like lfl9ll . 
The rates are spatially invariant except for a small region near 
the edge of the condensate band where the scattering phase 
space in the harmonic trap is modified by the energy cutoflF. 
Since the variation with radius is usually not significant, for 
our simulations we will use the approximate growth rate 



G(x) a! 7 s G'l" + G^' 



(All) 



which is exact in the inner region and a good approximation 
elsewhere. 



2. Scattering 

Our starting point of the scattering amplitude calculation, 
that is related to the actual scattering via Fourier transform, is 
Eq. (86) ofRef. iH 



>Rnc \^'^> -JRnc 

xF(x,Ki)[l -hF(x,K2)], 



(A12) 



where Ai2(e, k) = 6(u)i + u)2 - e/S)5(Ki + K2 - k). Upon 
evaluating the momentum delta-function we find 

M(x,k,0)=— — «f3KF(x,K)[l+F(x,K)] 

x5{2 |k - ^fixxj -k-k^j. (A13) 
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FIG. 12: Spatial variation of Gi{Q) in a spherical trap V(Q) = 
mLtr^Q~/2. Growth rates are shown for yu = O.lkgT (solid line) and 
yu = O.OlksT (chain line). The dashed vertical line indicates the point 
V(Q) = 2E]t/3 where the rates begin to vary. The condensate band 
terminates at 2^ = '\j2Eii/nnoj. The ratio Rq (see text) is shown 
versus j^/k^T (inset). 



Using the notation Ha)(x, K) = eK-mSixx//!(x), where 

ckCx) = + Veff(x), (A14) 

2m 

in the variable q = K - m(Q x x)/H we have 

M(x,k,0) = —-3 A<5 2q-k-k2 , 

(2;r)';z-' Jr„^ ^ ^ 

x/r,^(eq(x))[l + /r,;.(eq(x))], (A15) 

where fT,^,(z) = [e'---t'^"">^ - Writing M(x,k,0) = 

Ml (x, k)+M2(x, k) to distinguish the term with one or two BE- 
distributions as before, and choosing the z-component of Ki 
along k, we find, for the term involving one BE-distribution 



Mi(x,k) = 



4mu^ n 
(2WW ^ 

X I KxdKie-'''^^''^y-"'®mi2Kx < 1). 
JKk(:. 



(A16) 



The theta function generates the effective cutoff Ki > 
Kmini^) = ^3x{Kr{x), |k|/2). Evaluating the integral gives 



Mi(x,k) = 



nm 

{2nffP- j6ft2|k| 



(A17) 



2' 



which can be further simplified as follows. We want the 
condition for K,„j„{x) - Kr{x), which can be expressed as 
^^k-/8m < Er-Vss{x). The restriction of |k| is determined by 
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the requirement that k = ki + k2, where fp^k^^/lm + Veif (x) < 
Er for atoms in the condensate band. Then ffh?/^m = 
n\\/^m + n\l/^m + n\i ■k2/4m < - VeffCx) and the con- 
dition always holds. The exponential argument in Eq. ( IA17b 
then simplifies to /3{Er - jj). The term involving two BE- 
distributions, M2(x, k), is computed similarly, and finally we 
arrive at 



M(x, k, 0) = 



M 



(27r)3;!|k| (e/S(£«-//) _ 1)2 " (27r)3|k| 



(A18) 



As the result is independent of x, in this paper we drop the 
redundant arguments and in place of M(u, v, 0), we obtain 
the the scattering amplitude given by the Fourier transform 
ofEq. (IaTsT i 



M(v) 



M r 

(2n)^ J 



-iky 



(A19) 



APPENDIX B: NUMERICAL ALGORITHM FOR THE PGPE 
IN A ROTATING FRAME 

In this appendix we outline our numerical method for 
evolving the two-dimensional SGPE in cylindrical coordi- 
nates in a rotating frame of reference. For our purposes 
the method has several advantages over Fourier or Crank- 
Nicholson methods. It is a pseudo-spectral method makes use 
of a mixture of Gauss-Laguerre and Fourier quadratures for 
efficiency, and most importantly, allows an energy cutoff to be 
efficiently implemented in the rotating frame. 

Since the noise in the simple growth SGPE is additive, we 
use the stochastic Runge-Kutta algorithm due to the superior 
stability properties of Runga-Kutta methods ll46ll . Conse- 
quently, our central requirement is a means to efficiently eval- 
uate PLcpaix, r), or the projected Gross-Pitaevskii equation 



da 



(Bl) 



where is the two dimensional Laplacian. The determinis- 
tic part of the SGPE evolution is then found by transforming 
to an appropriate rotating frame. The solutions of the linear 
Schrodinger equation 



1. 



r2 



- Vi + indg + - Y„,(r, 6) = CJ^,Y„i{r, 0), 



(B2) 



are the Laguerre-Gaussian modes 



Yni(r,0) 



\nin + \l\)\ 



(B3) 



where L"'(x) are the associated Laguerre polynomials, and n 
and / are the quantum numbers for radial and angular excita- 
tions. The spectrum of the modes is 



wf, = 2n + |/| -0/+ 1, 



(B4) 



where « e {0, 1, 2, . . . } / e {. . . , -1, 0, 1, 2, . . . ). 

Imposing an energy cut-off in the rotating frame causes the 
range of available energy states to be biased toward states ro- 
tating in the same direction as the frame. We define the energy 
cut-off in the rotating frame by keeping all modes satisfying 



2n- 



■D.l<N, 



(B5) 



corresponding to cut-off energy Er = tiUriN + 1). From Eq. 
\, the quantum numbers are restricted by the cutoff to 



< « < [a^/2] 
-L{n) < I < h(n). 



(B6) 
(B7) 



where in this appendix [x] represents the integer part of x, and 
l±{n) - \{N - 2n)(l + Q)j. The projected wave-function then 
takes the form 



[N/2] ;+(«) 

«=0 l=-Hn) 



(B8) 



Projecting ( IBll i onto the single particle basis (|B3t leads to the 
equation of motion 



I—— = a>„,a„i + AF„i(a), 

OT 



(B9) 



where 

Fni{a) 



r-2n /-»oo 

= I de \ 

Jo Jo 



rdr Y*i(r, e)\a(r, 0)\^a(r, 0) (BIO) 



is the projection of the GPE nonlinear term, which is the nu- 
merically difficult term to evaluate for any pseudo-spectral 
method of the solving the GPE. Changing variables gives 



2 r^ZK y-»O0 

Fni(a) = — dOe-"" dx^Jx)7r 
27!- Jo Jo 

X\a{ 0)\^a{ V^, 6*) 



(BID 



where ^niix) = Y„i{ y/x, 0)e The integral to evaluate at 
every time step of propagation is always of the form 



1 p2n y-»D( 

- d0e-"' 
2^ Jo Jo 



dxe-^''Q{x,0) (B12) 



where Q{x, 0) is a polynomial in x and e'^ of order determined 
by the cutoff. The angular integration can therefore be car- 
ried out efficiently using an EFT. The order of e'" in a(x) is 
constrained by ^} to be -L = -/_(0) < / < /^(O) = l+. 
Choosing the grid so as to accurately compute all integrals 
^ d0 e^'''', where q is an integer in the range -2(L -i- 1+) < 
q < 2(L + l+), the integral can be carried out exactly for the 
condensate band by discretising as 0j - jA0 = jln/Ng for 
= 0, 1, . . . , A^fi - 1, where Ng = 2(L + h) + 1, to give 

d0 e-"'Q(x, 0) = X dj). (B13) 

j=0 
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The X integral is computed by noting that the polynomial part 
of the wave-function is of maximum order [^/(l - Q)j /2 in 
X (since x^'''l~L!j}{x) is a polynomial of degree n + |/|/2 < |/+|/2) 
and so the integrand is of maximum order 2|^^/(1 - Q)j , 
which can be computed exactly using a Gauss-Laguerre 
quadrature rule of order Nx = \n l{\ - Q)j + 1. The radial 
integral in (IB 12b is 



dx Q(x, e)e 



-2x 



1 ^' 

S 

k=\ 



WkQ(xkl2,e), (B14) 



where jc/. are the roots of Lff^{x) and the weights are known 



constants 114711 . We can thus write (IBl lb as 



F,Aa) = y m <^ni(xk/2) —Ye 



k=l 



(B15) 



x|ff( y/^,0jra{ y/^,6j). 



where the combined weight variable is % - 7rwAe '"'72. Once 
the modes are pre-computed at the Gauss points 



kn 



<J>„\i\(xk/2) 



(B16) 



the mixed quadrature spectral Galerkin method for evolving 
the GPE is 

(i) Compute the radial transform 

lN/2] 

Xkl - Yj ^^kn ^nl- (B17) 
11=0 



(ii) Construct the position field at the quadrature points 
using the FFT, after zero padding Xki to length Nq in the 
/ index, jpi/ = pad,(;t'«> A^e): 



kj 



(B18) 



(iii) Compute the nonlinear term 

^kj = \'Vkj?^kj. 



(B19) 



(iv) Compute the IFFT 



N„-\ 



®- = We 2 



(B20) 



(v) Calculate the Gauss-Laguerre quadrature 



F„iia) = J] m ®u pI,. 



(B21) 



We emphasize that the transform matrix f [^.J^ is separable, 
but requires a different kyi matrix for each |Z|. The extra mem- 
ory required to store (cf. the method of ll48ll for the non- 
rotating GPE) is mitigated by the modest memory used by the 
FFT. The order of the FFT rule required is consistent with the 
one dimensional harmonic oscillator rule given in [48;] where 
it was shown that to numerically evolve a Gaussian weighted 
polynomial wave-function of order .^V in x according to the 
PGPE, a Gaussian quadrature rule of order 2N + 1 is required. 
Since the FFT is a quadrature rule for unweighted polynomi- 
als in e'" we can expect the same behavior. This is evident 
once we note that L + 1+ is the order of the wave-function in 
powers of e'" and therefore plays the same role as in the 
standard Gauss-Hermite quadrature. 

The discretisation of angle and radius have thus been cho- 
sen as the smallest number of points needed to faithfully 
project the nonlinear term onto the representation basis with- 
out numerical error. The computational load of the method 
scales as O(NNxNg) where = max(A^^-, log Ng), affording 
a significant advantage over Cartesian grid methods in the 
rapidly rotating regime since Cartesian grids are compara- 
tively inefficient at representing the quantum states involved. 
Most notably, in the limit Q — » we may truncate the basis 
at the lowest Landau level Nx - 1 and the method reaches an 
optimal efficiency of one dimensional FFT. 
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